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ABSTRACT 

We investigate the evolution of intracluster medium during a cluster merger, 
explicitly considering the relaxation process between the ions and electrons by 
N-body and hydrodynamical simulations. When two subclusters collide each 
other, a bow shock is formed between the centers of two substructures and 
propagate in both directions along the collision axis. The shock primarily heats 
the ions because the kinetic energy of an ion entering the shock is larger than 
that of an electron by the ratio of masses. In the post-shock region the energy is 
transported from the ions to electrons via Coulomb coupling. However, since the 
energy exchange timescale depends both on the gas density and temperature, 
distribution of electron temperature becomes more complex than that of the 
plasma mean temperature, especially in the expanding phase. After the collision 
of two subclusters, gas outflow occurs not only along the collision axis but 
also in its perpendicular direction. The gas which is originally located in the 
central part of the subclusters moves both in the parallel and perpendicular 
directions. Since the equilibrium timescale of the gas along these directions 
is relatively short, temperature difference between ions and electrons is larger 
in the directions tilted by the angles of ±45° with respect to the collision 
axis. The electron temperature could be significantly lower that the plasma 
mean temperature by ~ 50% at most. The significance of our results in the 
interpretation of X-ray observations is briefly discussed. 

Subject headings: galaxies: clusters: general — hydrodynamics — intergalactic 
medium — plasmas — X-rays: galaxies 

1. INTRODUCTION 

Clusters of galaxies (CG) contain collisionless particles, galaxies and dark matter, and 
diffuse gas component. The gas component is called intracluster medium (ICM). The 
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ICM is the plasma with temperature of about 10 8 K, thus emitting X-ray mainly through 
the thermal bremsstrahlung of the electrons (Sarazin 1988). Since the ICM is optically thin 
for X-rays and spreads over the entire cluster, it can provide good information regarding 
the physical property of the clusters. 

According to the hierarchical clustering scenario, it is believed that CG are formed 
through subcluster mergers and/or absorption of smaller galaxy groups. The rate of mergers 
are dependent on the cosmological parameters, especially on the density parameter, Q 
(Richstone, Loeb, & Turner 1992). However, when this is applied to more realistic situations 
such as observational results (eg., West & Bothun 1990; Rhee et al. 1991; Jones & Forman 
1992; Mohr, Fabricant, & Geller 1993; Bird 1994; Serna & Gerbal 1996; Gurzadyan & 
Mazure 1998; Solanes, Salvador-Solo, & Gonzalez-Casado 1998) or cosmological numerical 
simulations (Evrard et al. 1993; Jing et al. 1995; Crone, Evrard, & Richstone, 1996), 
uncertainty of dynamical property of CG during mergers seriously affects the estimation of 
Q (Nakamura, Hattori, & Mineshige 1995). 

To study evolution of CG, plenty of numerical simulations have been performed so 
far using N-body and hydrodynamical codes. Perrenod (1978) was the first to calculate 
the evolution of a spherical symmetric CG with a standard mesh hydrodynamical code. 
Recently, Takizawa& Mineshige (1998) calculated the spherical CG and investigate the 
effect of pre-heating on ICM density profiles. To investigate the formation and evolution in 
more realistic situations, many simulations starting from the cosmological initial density 
perturbations are performed (e.g., Evrard 1990; Katz & White 1993; Metzler & Evrard 
1994; Navarro, Frenk, & White 1995; Bryan & Norman 1998; Eke, Navarro, & Frenk 
1998; Yoshikawa, Itoh, & Suto 1998). Schindler & Miiller (1993) find that characteristic 
temperature structures occur in merging clusters through shock heating and adiabatic 
compression and expansion. On the other hand, there are simulations starting from rather 
idealized initial conditions to study the mergers in detail. Roettiger, Loken, & Burns (1997) 
study general behaviors of head on collisions with a high resolution finite-difference code. 
Their observational consequences are discussed in more detail in Roettiger, Burns, & Loken 
(1996). Ishizaka (1997) find that the specific energy ratio of ICM and galaxies, (3 spec , is a 
good indicator of merging clusters and that it can be used to determine the phase of the 
merger. In particular, off-center collisions are investigated by Ricker (1998). 

In recent years, rather complex electron temperature structures are revealed in merging 
CG through X-ray observations (e.g., Fujita et al. 1996; Honda et al. 1996; Churazov et 
al. 1998; Donnelly et al. 1998a; Donnelly et al. 1998b; Davis & White 1998; Markevitch, 
Sarazin, & Vikhlinin 1998). In addition, Hanami et al. (1998) find that the energy of 
iron K a line is higher than that expected from the temperature in ionization equilibrium 
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plasma in SC1329-313 in Shapley concentration. These results implicate there exits strong 
bulk flow motion in ICM of merging clusters. Temperature maps of ICM contain important 
information of mergers as well as f3 spec does. Some numerical simulations are tried to 
explain these observational results; Roettiger, Burns, & Pinkney (1995) tried the modeling 
for A2256, and Ishizaka & Mineshige (1996) for Coma. More detailed "anatomy" is done 
for A754 by Roettiger, Stone, & Mushotzky (1998). 

As mentioned above, temperature structures of ICM provide us with very important 
information about the phase of the mergers. However, we can find only the electron 
temperature through X-ray observations until now. In all the previous studies the discussion 
is based on the assumption that observed electron temperature (T e ) is the same as the 
plasma mean temperature (T). However, we believe this assumption is problematic. The 
shock primarily heats the ions because the kinetic energy of an ion entering the shock is 
larger than that of an electron by the ratio of masses. In addition, the equilibrium timescale 
between electrons and ions (i e q) in typical CG becomes 10 8 ~ 10 9 yr if Coulomb coupling is 
considered as the relaxation process (see Spitzer 1962): 

t eq = 2.0 x 10 8 yr -! —?-) , 1 

q 3 V 40 J V10- 3 cm- 3 / V10 8 K7 ' v ; 

where n; is the ion density, T e is the electron temperature, and In A is the Coulomb 
logarithm. Indeed, some CG have larger electron temperature gradients than those expected 
in numerical simulations (Markevitch 1996). A two-temperature model of ICM is one of 
solutions to solve this discrepancy (Fox & Loeb 1997; Chieze, Alimi, & Teyssier, 1998; 
Ettori&Fabian 1998; Takizawa 1998). In cluster mergers, it is believed that shock heating 
plays an important role in the evolution of ICM. Furthermore, according to the previous 
studies about mergers, timescale surviving complex temperature structures is an order 
of 10 9 yr, which is comparable to t eq . Therefore, difference between T e and T should be 
considered properly to study the evolution of T e structure in mergers. 

To investigate the evolution of two-temperature ICM in cluster mergers, we perform 
N-body and hydrodynamical simulations considering the relaxation process between 
ions and electrons. For the hydrodynamical part, we choose the smoothed particle 
hydrodynamics (SPH) method. Although SPH codes are not better for problems with 
shocks than mesh codes, its Lagrangean nature is suitable to deal the electron-ion coupling 
simply. 

The rest of this paper is organized as follows. In §Q we estimate analytically some 
physical quantities relevant to cluster mergers. In §| we describe the adopted numerical 
methods and initial conditions for our simulations. In §^ we present the results. In §[| we 
summarize the results and discuss their implications. 
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2. BASIC CONSIDERATIONS OF EQUILIBRIUM TIMES 

In this section, we consider the scaling laws of CG and the initial conditions of cluster 
mergers. We assume that the power spectrum of the cosmological density fluctuation field, 
P(k), obeys power law, P(k) oc k n . According to the spherical model for the non-linear 
collapse (see Peebles 1980, Padmanabhan 1993, etc.), the virial radius (r), the density (p), 
and the virial temperature (T) of the dark halo obey the scaling laws as follows, 

r oc M (n+5)/6 , (2) 
pocM-(" +3 )/ 2 , (3) 
T oc M (1 " n)/6 , (4) 

where M is the total mass of the halo. The effective index of P(k) of CDM is —1 ~ —2 in 
the galactic and CG scale. We consider the case of n = —2. Then the above relations are, 

r oc M 1/2 , (5) 
p oc M- 1/2 , (6) 
T oc M 1/2 . (7) 

Using these relations and the gas fraction, f g , we obtain the scaling laws of equilibrium 
timescale between ions and electrons, t cq oc (/ g p) _1 T 3 / 2 , as follows. 

t eq oc f-'M^. (8) 

Thus, the more massive CG is, the longer t eq is, if f g is constant. However, the observational 
results about the luminosity-temperature (LT) relation suggests f g depends on M 
(Mushotzky 1984; Edge & Stewart 1991; David et al. 1993; Fukazawa 1997; Mushotzky & 
Scharf 1997; Markevitch 1998; Arnaud & Evrard 1998). Since L oc (/ g p) 2 T 1 / 2 J R 3 , / g oc M 3 / 8 
if we assume L oc T 3 . Then we obtain 

t eq oc M 7 / 8 . (9) 

t eq is longer in more massive CG also in this case, though the dependence is slightly weaker 
than the case when f 9 is constant. 

Next, let us consider the the merger of subclusters whose masses are M x and M 2 , 
respectively. The collision velocity V is estimated as follows. From the dynamical energy 
conservation law, we obtain, 

GMM GMM l w . 



2R r 1 + r 2 2 
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where G is the gravitational constant, r\ and r 2 is the virial radius of the each subcluster, 
R is the maximum expansion radius of the whole system, and M = M 1 M 2 /(M 1 + M 2 ) is 
the reduced mass. Since the maximum expansion radius is twice the virial radius, using the 
relation (^) we obtain 

-R = 2x — — — n. (ll) 



V M 1 

It is convenient to introduce the parameter a = M 2 /Mi. Without losing the generality, we 
can set < a < 1. Then, from the scaling relations of @ and equations ( fLOD and ([11]), we 
can describe V as follows, 



l 



.1 + «V 2 4(l + a) 1 /2 
On the other hand, the sound speed of subcluster 1, c\ is 



(12) 



v GMi 

^ = 7(7-1)^, (13) 



where 7 is the ratio of specific heat. Thus, the Mach number, M.\ = V/ci, is, 



7(7 - 1) 



(14) 



.l + aV2 4(l + a)V2_ 
From the scaling relation (^|), the Mach number for subcluster 2, M.2-, is, 

A<2 = a -l/2 M 2 v (15) 

Figure || shows M.\ (solid line) and 7Vf 2 (dotted line) with respect to a when 7 = 5/3. 
For major mergers (a > 0.1), these values are greater than ~ 2. Thus, in CDM universe, 
supersonic collisions are quite natural in major mergers. 

The Mach number gives the maximum temperature difference between ions and 
electrons in the post-shock regions. The specific inertial energy of ions and electrons are, in 
the limit of no mutual energy exchange, 

^ ~ c 2 + V 2 , (16) 
fj,m p 

(17) 



ksT e 2 



where \i is the mean molecular weight, m p is the proton mass, k-Q is the Boltzmann constant, 
and c is the sound speed in the pre-shock regions. Thus the difference between T e and 
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f ~ (Ti + T )/2, is, 



T e 1 



(18) 



T l + .M 2 /2' 



(19) 



Thus, from figure [I], we expect T e /T is less than about 0.6. We will confirm this later by 
numerical simulations. 



3. 



THE SIMULATIONS 



In the present study, we consider the CG consisting of two components: collisionless 
particles corresponding to galaxies and DM and two-temperature gas corresponding to 
the ICM. When calculating gravity, both components are considered, although the former 
dominates over the latter. Radiative cooling and heat conduction are not included in our 
simulations. 



To solve the hydrodynamical equations for the gas component, we used the smoothed- 
particle hydrodynamics (SPH) method (see Monaghan 1992). Although more accurate 
finite-difference methods better handle shocks than SPH methods, its fully Lagrangean 
nature is convenient to calculate the evolution of two-temperature plasmas since the 
coupling between ions and electrons can be treated more simply in Lagrangean view as 
follows. 

As the standard SPH codes for one-temperature fluid, we solve the continuity equation, 
the momentum equation, and the thermal energy equation with artificial viscosity to treat 
the shocks. In addition to these equations, we solve one more equation for the normalized 
electron temperature, T e = T e /T, where T e is the electron temperature and T is the plasma 
mean temperature. We assume that artificial viscous heating is effective only for ions 
and that only the Coulomb coupling is considered in the relaxation process. Then, the 
Lagrangean time evolution of T e is (see Appendix), 



where 2] = Xj/T is the normalized ion temperature, Q V { S is the artificial viscous heating 
per unit mass, and u is the thermal energy per unit mass. For the numerical integration 



3.1. 



The Numerical Method 




(20) 
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of equation (p0|), we should take a time-step shorter enough than not only the Courant 
and viscous timescale as usual SPH codes (see Monaghan 1992 §10.3) but also t eq . In our 
simulations, however, the latter is much shorter than the former in the central high-density 
regions. Thus, the integration with timestep relevant to t eq takes a huge computational 
time. In addition, since two-temperature nature is not important in such regions, the 
time-step control above mentioned is not a good choice. To integrate equation (|20|) with 
the Courant and the viscous time-step control, we use the results of Fox & Loeb (1997) 
as follows. Let us integrate equation ( ^20]) from t = to to t = to + At. First, we integrate 



equation ( p0|) regarding the second term on the right-hand side being negligibly small. In 
this case, neglecting the small change of Coulomb logarithm in t eq , we can integrate this 
analytically (Fox & Loeb 1997), 



At m 



t 2s Tli + Ue . 

where, 



F(f c (t + At))-F(f e (t )) 



(21) 



^=Ki^)-^( 1+ !)' ^ 



and t 2s = t eq (to)T e (to)~ 3 ^ 2 . Second, we add the contribution from the second term to 
T e (to + At) using the second-order Runge-Kutta method. If the time step is chosen correctly, 
the displacement of the position of each SPH particle is smaller enough than the spatial 
resolution and Q v i S /u is smaller enough than unity. Thus the above method can follow the 
evolution of T e efficiently and with reasonable accuracy. 

Gravitational forces are calculated by the Barnes-Hut tree algorithm (Barnes & Hut 
1986) and softened using a Plummer potential profile $ oc (r 2 + e 2 )" 1 , where e is the 
softening parameter. We set e one tenth of the initial core radius of the smaller subcluster 
in the simulation. Tree structure is also used to search for nearest neighbors in SPH 
calculations (Hernquist & Katz 1989). In the standard one- dimensional shock tube test, 
our code nicely reproduced the analytic solution. The shock and contact discontinuity is 
broadened over a range of about three times of smoothing length. During the calculations 
described below, the total energy is conserved within ~ 0.4%/(1000 steps). 



3.2. Models and Initial Conditions 



We consider mergers of two virialized subclusters of galaxies with masses M\ and M 2 . 
Initial configuration of each subcluster is described as follows. The spatial distribution of 
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the DM in each subcluster is represented by the King distribution; 

PDM(r) = \0 (r>r out ) (23) 

where, po is the central density, r c is the core radius, and r out is the radius of the cluster. 
We set r ou t = 5 x r c . This is nearly equal to the virial radius of the King distribution. p is 
determined by the total mass of DM (M DM ), 

U ■ 1" 1 



== + hi|t ou t + Vw + X l 



(24) 



where t out = r out /r c . We assume velocity distribution of the DM particles to be an isotropic 
Maxwellian. Then, from the second moment of the collisionless Boltzmann equation, 
one-dimensional velocity dispersion at r, (i\{r), is (e.g., Binney & Tremaine 1987) 

2/ x 1 r° ut / x GM r , , v 

0i r = TT Pdm r — (lr, 25 

Pdm(h Jr r 2 

where M r is the mass inside r. We assume that initial ICM temperature is isothermal and 
equal to the virial temperature at r out . The ICM is initially in hydrostatic equilibrium 
within a cluster potential of the DM and the ICM, itself. 

We set the initial conditions as follows. Two subclusters are placed at rest with a 
separation between the center of each subcluster, r sep , initially. The coordinate system is 
taken in such a way that the center of masses is at rest in the origin. Two representative 
cases are examined in the preset study: a collision between two equal-mass clusters (Run 
A), an absorption of a smaller cluster by a larger one (Run B). The scaling law between 
the two subclusters and r sep are determined as the case of P(k) oc k~ 2 in In Run 
A, two subcluster have the same masses; Mi = M 2 = 0.5 x 1O 15 M . We set the core 
radii to be r c = 0.2Mpc. Each subcluster consists of 5000 collisionless particles and 5000 
SPH particles. The total gas mass fraction is 10% for both clusters. In Run B, the mass 
ratio is M\ : M 2 = 4:1. The larger cluster has the same mass as that in Run A. The 
particle numbers of the smaller one are one fourth of the larger one. The parameters in our 
calculations are summarized in table |l|. 



4. RESULTS 

4.1. Time evolution of the total feature during the mergers 



Figure |2| shows the time evolution of various physical quantities representing the 
total feature of the system of Run A: from the top, the total luminosity (L), the 
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emissivity-weighted mean electron temperature (T e ew ), the square root of the velocity 
dispersion parallel to the collision axis {&//), and the specific energy ratio of DM and ICM 
(Ajpec)- Pspec is defined as, 

&pec = ^f^- (26) 

When calculating L, we assume the thermal bremsstrahlung of optically thin plasma 
(Rybicki & Lightman 1979) and neglect the line emissions. 

Around the collision, all of L, T 6jew and an rise but on different timescales. Since the 
rise of T e ew is quicker than that of a//, /3 spec has two peaks in the time evolution. The 
minimum of /3 spec is coincident with the maximum of T e ew . 

Similar results are seen also in Run B (figure [7]). However, the amplitudes of time 
fluctuation are smaller than in Run A. 

Since the emissivity is proportional to the square of the gas density, emissivity-weighted 
mean temperature (EWMT) is more like the temperature in the central high density regions, 
where t eq is ~ 10 8 yr or less. Thus, the difference in EWMTs between electrons and ions 
is practically negligible. Indeed, this difference becomes ~ 1% even at maximum in both 
runs. Therefore, our results agree qualitatively with those of the previous one-temperature 
simulations such as Ishizaka (1997). 



4.2. Time evolution of temperature distribution of Run A 



As seen in j j4.1| , two-temperature nature of ICM does not manifest itself in EWMT of 
CG. This is not the case, however, when we can resolve the spatial distribution of ICM 
temperature. Let us examine first the case of Run A. 

First of all, let us see the T map and the gas velocity field, and discuss the dynamical 
evolution of Run A. Figure |3| shows the snapshots of X-ray surface brightness (contours) 
and emissivity-weighted T (colors) distribution seen from the direction perpendicular to 
the collision axis. X-ray surface brightness contours are equally spaced on a logarithmic 
scale and separated by a factor of 7.4. The blue, green, yellow, and red colors correspond 
to k?,T ~ 5 keV, 10 keV, 15 keV, and 20 keV, respectively. The times are listed above 
each panel: t = 4.1, 4.3, 4.6, and 4.85 Gyr correspond to the first maximum of /3 spec , the 
minimum of /3 spec , the second maximum of f3 spec , and the minimum of T 6jew , respectively. 
Figure ^ shows the evolution of the gas velocity field, seen from the same direction at the 
same times as in figure |^. The longest vector corresponds to the maximum velocity listed 
below each panel. 
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When two subclusters just contact each other at t — 3.5 Gyr (upper left panel), T 
rises slightly at the interface between the two subclusters due to the adiabatic compression 
but we can clearly distinguish two individual clusters through the X-ray image. Then the 
two approach each other and we see double peaks in the X-ray emissivity profile of 'one 
cluster' at t — 4.1 Gyr (upper middle). Just between the peaks T rises up to ~ 20 keV 
due to a shock. The shock front is perpendicular to the collision axis. However, since 
T still remains nearly the initial value (~ 5 keV) around the peaks of the X-ray image, 
EWMT does not increase so much that /3 spec becomes rather large (figure §). At the most 
contracting phase (t = 4.3 Gyr: upper right), two peaks merge to one peak in the X-ray 
image and the X-ray image elongates to the directions perpendicular to the collision axis. 
Since the high temperature (~ 20 keV) region is located around the X-ray peak, EWMT 
becomes maximum (figure |2|). Then the cluster expands and two shocks propagate in the 
opposite directions along the collision axis (t = 4.6 Gyr: lower left). We can clearly see two 
'lens-shaped' high temperature regions associated to the shocks in the T structure. We can 
also see the cooler region (~ 5 keV) spread in the directions perpendicular to the collision 
axis. These characteristic temperature structures reflect the fact that the gas outflow occurs 
not only along the collision axis but also in its perpendicular directions (see figure ||). In 
these directions the gas is effectively cooled by the adiabatic expansion and emits more 
X-ray than the hottest regions. Thus, EWMT decreases, though high temperature regions 
still exist. At the most expanding phase (t = 4.85 Gyr: lower middle), the X-ray image 
elongates along the collision axis. The shocks reach in very low density region. We see 
an 'X-shaped' region with slightly high temperature (indicated by light blue color). This 
structure is also due to adiabatic expansion associated to the gas outflows. Then the system 
contracts again and settles down to the spherical structure (t = 7 Gyr: lower right). Owing 
to the subsonic accretion of the gas along the collision axis, slightly high temperature 
regions are created on the both sides of the collision axis (see figure |]). 

Note that there can be complex temperature structures in ICM when X-ray images 
do not have definite substructures. It is certain that image elongation is occurred in these 
case, but this is severely affected by the viewing angles. Thus, temperature map is more 
suitable for the diagnostic of merging clusters. 

Next, let us see the T e map and discuss the difference between T and T e . Figure [| 
shows the same as figure |3|, but for T e . Although in the contracting phase (t < 4.3 Gyr) T e 
maps are similar to the T maps, in the expanding phase (t > 4.3 Gyr) both maps are rather 
different. To see the difference between T and T e in detail, we make the snapshots of the 
distribution of T e = T e /T. Figure | shows the T e maps of Run A. The red, yellow, green, 
and blue colors correspond to T e ~ 0.1, 0.3, 0.5, and 0.7, respectively. 
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In the contracting phase (t < 4.3 Gyr), the difference between T and T e is rather small. 
At t = 4.1 Gyr, T e is slightly lower than T at the shock. In overall there are no significant 
two-temperature regions. Then, at t = 4.3 Gyr, split two-temperature regions appear due 
to the propagation of shocks. Since t eq is longer in the outer parts where density is low, 
two-temperature regions are more spread there than in the central parts. In this phase, 
however, two-temperature nature is less important than in the later phase because the the 
shocks is nearly standing shock and t cq is shorter in the central parts. 

On the other hand, in the expansion phase (4.3 Gyr < t < 4.85 Gyr), the distributions 
of T and T e are clearly qualitatively different. From figure |6] we find two common features 
both at t — 4.6 Gyr and at t = 4.85 Gyr. One is that T e is significantly lower (~ 50 %) than 
T in the post-shock high T regions. These regions are spread on ~ 0.5 Mpc scale behind 
the shocks. The other is that the difference between T and T e is larger in the directions 
tilted by the angles of ±45° with respect to the collision axis. This is also due to the gas 
outflows both along the collision axis and in its perpendicular directions. The gas originally 
located in the central part, whose t eq is relatively short, moves outward in these directions. 
Note that t eq is hardly changed by the adiabatic compression and expansion when 7 = 5/3. 
Thus, the difference between T and T e is less in these directions. As a result, we see four 
peaks in the T e map at t = 4.6 Gyr and the 'X-shaped' high temperature region becomes 
less contrasted in the T e map at t — 4.85 Gyr. Note that high temperature regions around 
~ 20 keV cannot be seen in the T e maps of the expanding phase. Therefore, it is probable 
that we underestimate the collision velocity by ~ 30% if we estimate it from the hottest 
regions in T e map in the expanding phase. Two-temperature nature still remains at t — 7 
Gyr in the outer parts as the case of spherical models (Fox & Loeb 1997; Takizawa 1998). 



4.3. Time evolution of temperature distribution of Run B 

First, we describe the dynamical evolution of Run B through the T map and the gas 
velocity field. Figure || shows the same as figure |3|, but for Run B. Temperature color scale 
is adjusted for Run B. The blue, green, yellow, and red colors correspond to k-^T ~ 4 keV, 
7 keV, 10 keV, and 13 keV, respectively. The times are listed above each panel: t = 3.9, 
4.1, 4.3, and 4.75 Gyr correspond to the first maximum of /3 spec , the minimum of /3 spec , the 
second maximum of /3 spec , and the minimum of T e ew , respectively. Figure |9] shows the same 
as figure f|, but for Run B. 

When the two subclusters approach each other, the bow shock with an arc shape is 
formed just between the two. Thus, the arc-shaped hot region is seen between the two 
peaks of the X-ray image at t — 3.9 Gyr (upper middle). This region clearly separates 
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the component of larger subcluster (~ 5 keV) from the component of smaller subcluster 
(~ 2.5 keV). Then, the two peaks merge to one triangle image of X-ray (t = 4.1 Gyr: 
upper right). The hottest temperature region (~ 13 keV) is located in the peak of the 
X-ray image and the hot region associated to the bow shock is seen elongated slantingly 
backward with respect to the motion of the smaller subclusters. Then the gas expands and 
the two shocks propagate outward along the collision axis. The one in front of the motion 
of the smaller subcluster (for the left of the figure) is arc-shaped and the other is rather 
flat (t = 4.3 Gyr: lower left). At the most expanding phase (t = 4.75 Gyr: lower middle) 
the X-ray image elongates along the collision axis and the shocks reach very low density 
regions. The elongation of the image is more significant for the left of the figure. Then the 
system contracts again and settles down to the spherical structure (t = 7 Gyr: lower right). 
Accretion flows are seen only from the front side with respect to the motion of the smaller 
subcluster. 

Next, we describe two-temperature nature of Run B through the T e map and the T e 
map. Figure shows the same as figure ||, but for T e and figure |TT| shows the same as 
figure H but for Run B. As the case of Run A, in the contracting phase (t < 4.1 Gyr) the 
T e distribution is similar to the T one though T e is slightly lower than T around the shock. 
On the other hand, in the expanding phase (t > 4.1 Gyr) substantial deviation emerges 
especially around the hot regions associated to the shocks. At t = 4.3 Gyr in the T map 
both two hot regions have almost the same temperature (~ 13 keV), whereas in the T e map, 
the former (with respect to the motion of the smaller subcluster) hot region becomes less 
than ~ 10 keV. On the other hand, in the backward hot region, ~ 13 keV component still 
remains. In this region, t eq is rather short because the gas there originates from the smaller 
subcluster (see the scaling relation ([8])). At t = 4.75 Gyr, as the case of Run A, T e is low in 
the directions tilted by the angles of ±45° with respect to the collision axis, although this 
tendency is less significant than in Run A. 



5. SUMMARY AND DISCUSSION 

We investigate evolution of ICM during mergers considering the relaxation process 
between the ions and electrons. From the simple analytical estimation, we find that in the 
CDM universe the equilibrium timescale between ions and electrons are longer in more 
massive CG even if the gas fraction is dependent on CG mass in such a way that the X-ray 
luminosity is proportional to the cube of the temperature. We estimate the collision velocity 
and show that supersonic collisions are quite natural in CDM universe in major mergers 
and that electron temperature can be less than a half of the plasma mean temperature in 
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the post-shock regions. Temperature estimated from X-ray observations can be significantly 
lower than ion temperature. 

We carry out the numerical simulations of the mergers by N-body and hydrodynamical 
simulations, incorporating the relaxation process between the ions and electrons. To solve 
the evolution of the normalized electron temperature, we adopt the results of Fox & Loeb 
(1997). The difference in emissivity- weighted mean temperatures of ions and electrons is 
practically negligible during the merger; the discrepancy is ~ 1% even at maximum. On 
the other hand, the spatial distribution of the electron temperature is significantly different 
from that of the plasma mean one especially in low-density regions in the expanding phase. 
In this phase, electron temperature is at most ~ 50% lower than the plasma mean one in 
the post-shock hot regions. In addition, the difference between them is more enhanced 
in the directions tilted by the angles of ±45° with respect to the collision axis. When 
the two subclusters have different masses, the hot region located to the former position 
(with respect to the direction of the motion of the smaller subcluster) has lower electron 
temperature than that located backward in the expanding phase. 

Recently, Markevitch et al. (1998) estimate the subcluster collision velocity of Cygnus 
A using the electron temperature map obtained by ASCA. From this temperature map and 
the ROSAT PSPC image, it is likely that Cygnus A is in just contracting phase such as Run 
A at t = 4.1 Gyr. Thus, two-temperature nature is almost negligible in this case. However, 
collision velocities are likely underestimated by ~ 30% if we use the electron temperature of 
the hottest regions of merging clusters in the expanding phase such as SC1329-313 (Hanami 
et al. 1998). 

Although we consider only head-on collisions for simplicity in the present calculations , 
off-center collisions should be investigated to model more realistic situations. Ricker (1998) 
find that that spiral bow shocks occur in off-center mergers, which probably produce rather 
complex electron temperature distribution. We will investigate this issue as a future work. 

We consider only the classical Coulomb coupling as the relaxation process between 
ions and electrons. It is possible, however, in ICM more efficient relaxation processes could 
be effective (McKee & Cowie 1977; Pistinner, Levinson, & Eichler 1996). In this case the 
equilibrium timescale could be shorter than the value given by equation (|I]). Therefore, the 
temperature difference between ions and electrons can be less than our results. If magnetic 
field exits in ICM, it is possible that electrons are also significantly heated by shocks 
associated with MHD instabilities. Also in this case the temperature difference between 
ions and electrons could be less. 

We neglect the heating process from the galaxies to ICM. If substantial amount of 
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thermalized hot gas is injected into ICM from the galaxies, temperature difference could be 
less. It is possible that mergers of CG activate the star formation in the member galaxies 
through the galaxy-galaxy interaction, galaxy-ICM interaction, etc (Caldwell et al. 1993) . 
In this case, this effect cannot be negligible but the effect on ICM temperature distribution 
is rather sensitive on the detailed modeling of the star formation activity (Fujita 1998). 

We adopt the SPH method for the hydrodynamical part in the present simulations 
for simplicity to treat the coupling between ions and electrons. However, resolution of the 
shock in SPH is not better than that in high resolution finite-difference codes such as TVD, 
PPM, CIP, etc, because artificial viscosity is used to treat the shocks in SPH. Furthermore, 
bulk viscosity is included in our code and possibly effects electron temperature slightly. 
Although we believe that our simulations can follow the qualitative behavior, it is useful to 
simulate two-temperature ICM during mergers with finite-difference codes. 

The author would like to thank C. Ishizaka, H. Hanami, and S. Mineshige for helpful 
comments. This work is supported in part by a Grant-in-Aid from the Ministry of 
Education, Science, Sports and Culture of Japan (6179) and Research Fellowships of the 
Japan Society for the Promotion of Science for Young Scientists. 



A. DERIVATION OF THE EQUATION FOR THE NORMALIZED 

ELECTRON TEMPERATURE 

In the SPH calculations shocks are treated through artificial viscosity. We assume 
that artificial viscous heating is effective only for ions and that only Coulomb coupling is 
considered in the relaxation process. Then the Laglangean time evolution of the electron 
temperature T e and the mean temperature T = (n e T e + n{T\)/{n e + n{) is 

dT e . ^.T e dn T[ — T e , . . 

if = < A1 > 

§ = + (A2) 

at n dt 3 k B 

where TJ is the ion temperature, n is the gas density, and 7 = 5/3 is the ratio of specific 
heats and Qvis is the artificial viscous heating per unit mass. We use the artificial viscosity 
in the form described in Monaghan (1992) §4.1, which is the most commonly used one. 
Then, the explicit expression of Q vis for the i-th SPH particle is 

Qvis,i = \ mj^ijivi - vj) ■ ViWij, (A3) 
1 3 
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where 11^ is given by 




—etc., 



for (vj 
for (vj 




r,-) < 0, 
r,-) > 0, 



Pij 



(A4) 



and 



hjjjvi -Vj)- (Vj - Vj) 



(A5) 



In the above expressions, rrij, Vj, and r*j are the mass, velocity, and position of the 
j-th SPH particle, respectively, c^, and /i^- are the average of sound speed, density, and 
smoothing length of particle % and j, respectively. We set a — 1, /3 — 2, and rf = 0.01, 
which are typical values in SPH calculations. VjWjj is the gradient of the j-th particle's 
kernel function at tv 

Introducing the temperature normalized by T, i.e., T c = (T c /T) and TJ = (Tj/T), we 



find 




(A6) 



where w = 3/cBT/(2/xm p ) is the thermal energy per unit mass. 
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Run A Run B 

M l /M 2 (10 15 M e ) 0.5/0.5 0.5/0.125 

r c ,i/r C)2 (Mpc) 0.2/0.2 0.2/0.1 

k^T x jk^T 2 (keV) 4.78/4.78 4.78/2.39 

r sep (Mpc) 4.0 3.3 

/ g 0.1 0.1 

e (Mpc) 0.02 0.01 

iVi/iV 2 (SPH) 5000/5000 5000/1250 

NJN 2 (DM) 5000/5000 5000/1250 



Table 1: Model parameters 
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Fig. 1. — The square of Mach numbers of two virialized clusters as functions of a (the ratio 
of the two cluster masses) during the merger for the case with P(k) oc k~ 2 . The solid and 
dotted lines show the values for the larger and the smaller cluster, respectively. 

Fig. 2. — The time evolution of various physical quantities representing the features of the 
entire system for Run A: from the top, (a) the total luminosity (L), (b) the emissivity- 
weighted mean electron temperature (T ejew ), (c) the square root of the velocity dispersion 
parallel to the collision axis [&//), and (d) the specific energy ratio of DM and ICM (/3 spec ). 

Fig. 3. — The snapshots of X-ray surface brightness (contours) overlaid with emissivity- 
weighted T (colors) distribution seen from the direction perpendicular to the collision axis 
for Run A. X-ray surface brightness contours are equally spaced on a logarithmic scale and 
separated by a factor of 7.4. The blue, green, yellow, and red colors correspond to k^T ~ 5 
keV, 10 keV, 15 keV, and 20 keV, respectively. The times are listed above each panel. 

Fig. 4. — The evolution of the gas velocity field for Run A seen from the same direction 
at the same times as in figure |||. The longest vector corresponds to the maximum velocity 
listed below each panel. 

Fig. 5. — Same as figure |], but for T e instead of T for Run A 

Fig. 6. — The snapshots of distribution of normalized electron temperature, T e = T c /T for 
Run A. The red, yellow, green, and blue colors correspond to T c ~ 0.1, 0.3, 0.5, and 0.7, 
respectively. 

Fig. 7. — Same as figure |2], but for Run B. 

Fig. 8. — Same as figure |3|, but for Run B. Temperature color scale is adjusted for Run B. 
The blue, green, yellow, and red colors correspond to k-^T ~ 4 keV, 7 keV, 10 keV, and 13 
keV, respectively. 

Fig. 9. — Same as figure f|, but for Run B. 
Fig. 10. — Same as figure |5], but for Run B. 
Fig. 11. — Same as figure H but for Run B. 
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